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ABSTRACT 

A phase-space generation algorithm, capable to efficiently integrate the squared 
amplitude of any scattering process, is presented. The algorithm has been im- 
plemented in a Monte Carlo program, PHEGAS, which, using HELAC, a helicity 
amplitude computational package, can be used for automatic cross-section 
computation and event generation. Results for several scattering processes 
with four, five and six particles in the final state are briefly presented. 



July 2000 



The study of multi-particle processes, like for instance four-fermion production in e + e~, 
requires efficient phase-space Monte Carlo generators. The reason is that the squared 
amplitude, being a complicated function of the kinemtaical variables, exhibits strong 
variations in specific regions and/or directions of the phase space, lowering in a substantial 
way the speed and the efficiency of the Monte Carlo integration. A well known way out 
of this problem relies on algorithms characterized by two main ingredients: 

1. The construction of appropriate mappings of the phase space parametrization in 
such a way that the main variation of the integrand can be described by a set of 
almost uncorrelated variables, and 

2. A self- adaptation procedure that reshapes the generated phase-space density in 
order to be as much as possible close to the integrand. 

Up to now such algorithms have been developed in several cases to deal with specific 
processes, like four-fermion jl], four-fermion plus a photon [pi [TI| |TTJ and six fermion |J 



production in e + e~ collisions, as well as in the framework of general-purpose computa- 
tional packages like CompHEP [|J and GRACE ||. It is the aim of this letter to present a 
generalized recursive algorithm, together with its implementation, that can be used for 
automatic cross-section computation for any multi-particle process. 

In order to construct appropriate mappings we note that the integrand, i.e. the squared 
amplitude, has a well-defined representation in terms of Feynman diagrams. It is therefore 
natural to associate to each Feynman diagram a phase-space mapping that parametrizes 
the leading variation coming from it. To be more specific the contribution of tree-order 
Feynman diagrams to the full amplitude can be factorized in terms of propagators, vertex 
factors and external wave functions. In general, the main source of variation comes from 
the propagator factors and therefore our aim is to construct a mapping that expresses the 
phase-space density in terms of the kinematical invariants that appear in these propagator 
factors. Since in principle we need as many mappings as Feynman diagrams for the process 
under consideration, we have to appropriately combine them in order to produce the global 
phase-space density. A simple and well studied solution to this problem was suggested 
some time ago in reference 0. Let us represent the normalized phase space-density of a 
mapping by a function where $ refers to the (3n — 4)-dimensional phase space, n 

being the number of produced particles. The overall density can be represented by 

M 
i=l 

where 

M 

0<Cti<l ^ OLi — 1 
i=l 
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and M is the total number of mappings. Since the result of the integration does not 
depend on the specific values of ctj, the so-called a priori weights, the latter can be used 
to optimize the Monte Carlo integration. A self-adaptation procedure therefore suggests 
itself: during the evaluation of the integral, aij are repeatedly redefined ||, so that the 
variance of the integrand is minimized. It should be mentioned however that other self- 
adapting approaches can be used as well 0. 

In order to describe the construction of the phase-space mappings, let us consider a 
typical process in which two incoming particles produce n outgoing ones. The phase space 
can be represented by 

i=l ) i=l 

where P = q± + q 2 with q%, q2 being the momenta of the incoming particles. 

A well known property of Eq.(|l|) is that the phase space can be decomposed as follows 

/ m dQ 2 \ 

d<S> n (P;p 1 ,p 2 ,...,p n ) = \Jl^Jd® m (P;Q 1 ,...,Q m ) 

) (2) 

where the subsets {r 1; r 2 , . . . , r ni } up to {si, s 2 , ■ ■ ■ , s Hm } represent an arbitrary partition 
of {pi,p2, . . . ,p n }. The above equation can be generalized recursively resulting in an 
arbitrary decomposition of d$ n . Feynman graphs can be seen as a realization of such a 
decomposition, this latter being identified with a sequence of vertices of the graph. For 
instance a three-particle vertex V — (Q — > Qi, Q 2 ) in a Feynman diagram can be seen as 
part of the phase-space decomposition 

d*„ = ...^^d* a (Q;Qi,Q a )... • (3) 

Z7T Z7T 

The appropriate sequence of vertices, {Vi,V 2 , . . . ,Vk} can be chosen in such a way that 
a recursive construction of the phase space is realized. For instance V\ should contain 
at least one incoming particle whose momentum is known. The rest of the sequence is 
chosen recursively: vertex Vj is characterized by an incoming momentum Q which has 
already been generated in one of the {Vi, . . . , V^-i} and outgoing momenta Qi and Q 2 
that are generated according to Eq.(|3]). 

As a more illustrative example let us consider the graph of Fig.[l] for the process 

e~(lx) e+ (<?2) -> H~(pi) Vv(P2) u(p 3 ) d(p 4 ) 
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The appropriate sequence of vertices can be chosen as 

Vi = (?i -92, Q), ^ 2 = (Q-> Qi,Qa), V" 3 = (Qi -^pi,p 2 ),^ 4 = (Qa ^P3,P4) 

where Q, Qi, Q2 are the momenta of Z, W 7- , W 74 " respectively. In the first vertex, V\, both 
gi and q2 are considered known so this is a mere definition of Q = q\ + qi- The rest of the 
sequence realizes the following phase-space decomposition 

d$4(<?i + <?2;Pi,P2,P3,P4) = d$ 2 (<?i + Q2] Qi, Q2) 

dQl dQ 2 2 

d$ 2 (Qi;pi,P2) 

^ 2 (Q 2 ;P3,P4) (4) 

allowing for the correct treatment of the Breit-Wigner propagators of in terms of the 
variables Q\ and Q\. 

In the general case one can distinguish two types of vertices: 

1. All outgoing momenta involved in the vertex are time-like. 

2. One of them is space-like. 

It is worthwhile to mention that for 2 —>■ n scattering these two cases are the only possible 
onesQ. 

1 Notice that in the present study we restrict ourselves to scattering processes whose amplitudes do 
not exhibit non-integrable singularities over the available phase space. 
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For the first case the phase space decomposition can be written as 



Q 




dQl dQ 2 2 
2?r 2tt 

dQj dQ\ 
2tt 2tt 



d*2(Q^Qi,Q 2 



dcosd 



X^iQ^QlQl 



32vr 2 Q 2 

X(x, y, z) = x 2 + y 2 + z 2 — 2xy — 2xz — 2yz 

with the understanding that whenever Qi or Q2 represents an external momentum the 
corresponding factor dQ 2 /27i is set to 1. Generation can now proceed straightforwardly, 
by first generating Ql and Ql using any prescribed density, as well as cos 8 and in the 
rest frame of Q. Then by using the known momentum Q a boost to the initial frame 
can be performed. As it is easily seen the first case results to a rather simple generation 
algorithm. 

The second case is more involved. The phase space is decomposed as follows: 



Q 



Q 1 -q 7 



with 



d$ n = ...^^d^iQ^Q^Qz)... 

= ...jgf jgl dtrf0 1 (5) 
2tt 2tt v 32tt 2 Q \q 2 \ V ; 

t = {Qi- g 2 ) 2 = m 2 2 + Ql- ^{Q 2 + Q\- Ql) + ^-\q 2 \ cos9 

and (E 2 , q 2 ) being the incoming momentum q 2 in the rest frame of Q. In order to have 
an efficient generation according to Eq.(0) we need to know the limits of the t— and 
Ql~ integration: a detailed presentation of their derivation can be found in the Appendix. 



5 



Although in the two cases described so far we used a three-particle vertex the algorithm 
can be generalized in a straightforward way in the case of a four-particle vertex, either 
using the three-body phase space explicitly 

^ = dQf ^| d^_ d ^ {Q ^ Q ^ 

Z7T Z7T Z7T 

in the case all momenta are time-like, or using 



d* n = ...^^d* i (Q-+Q 1 ,Q m ) 



followed by 

...^^ 2 (g 2 3-Q 2 ,Q 3 )... 

Z7T ATI 

in the case one space-like outgoing momentum is present. 

Following the above described algorithm we end up with an expression for the phase- 
space density, 

d$ n -> n n n # fe n <* cos ^ (6) 

where s, and £j refer to the kinematical invariants entering the propagator factors of the 
graph and <pk and cos#z represent center-of-mass angles needed to complete the phase 
space parametrization. It is now straightforward to generate s, and tj with a probability 
density given by: 

• p(x) ~ (x — m 2 ) 2 + m 2 T 2 for massive unstable particles, like W ± , Z, t . 

• p(x) ~ x v for time-like massless propagators, e.g. 7, gluons, 

massless fermions. 

• p{x) ~ \x\ v for space-like massless propagators. 

so that the corresponding propagator factor cancels out in the Monte Carlo weight. The 
value of the exponent i>, for 7 and gluons, is chosen very close to 1 in order to account 
for the leading single-pole behaviour of the squared amplitude as a result of the gauge 
cancellations. 

The implementation of this algorithm, called PHEGAS, is based on and combined with 
HELAC II a package that computes any tree-order matrix element. HELAC is based on the 
Dyson- Schwinger recursive equations that proved to be superior to the Feynman diagram 
representation for amplitude computation. On the other hand it is still an open problem 
how to use Dyson- Schwinger representation to define phase-space mappings. We have, 
therefore, implemented an algorithm that allows the construction of all Feynman diagrams 
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form the HELAC solution of the Dyson- Schwinger equations, in a from suitable to be used 
by PHEGAS. In fact each Feynman diagram is represented by a sequence of integer arrays 
corresponding to its vertices. The user supplies information concerning the process under 
investigation such as the flavours of incoming and outgoing particles as well as a couple 
of control parameters as described in reference ||, along with a user-prescribed routine 
that specifies the desired cuts on the kinematical variables. The HELAC-solution of the 
Dyson-Schwinger equations for the process under consideration is used to produce the 
appropriate representation of all Feynman graphs. This information is then introduced 
into PHEGAS which produces phase space points according to the parametrization suggested 
by the corresponding mapping as well as the appropriate weight, taking automatically 
into account the prescribed phase-space cuts. The global density is then constructed by 
computing phase-space densities for all mappings followed by a multichannel optimization. 
The output of the program provides the total cross section as well as any kinematical 
distribution prescribed by the user. 

In order to show explicitly the usefulness of the proposed algorithm we consider the 
following typical examples of cross-section computation. 

• e~e + — > {iv^ud 

This is a well studied process within four-fermion physics at LEP2. We present here 
results form PHEGAS/HELAC in comparison with results from EXCALIBUR || |TJ. They 
are summarized in the following table: 





MC points 
w > 


result 
(fb) 


error 
(fb) 


efficiency 

(%) ' 


PHEGAS/HELAC 


1510700 


608.64 


0.61 


3.5 


EXCALIBUR 


1574175 


608.22 


0.57 


3.6 



where we have used 2 x 10 6 MC points, at ^/s = 190 GeV, and a fixed width pre- 
scription for internal unstable-particle propagators in both programs and identical 
input parameters. In the last column of the table the efficiency of the generator 
is given. The efficiency of an event-generator is defined as the ratio of the mean 
to the maximum Monte Carlo weight and it is also related to the number of the 
unweighted events: for instance in the above run a sample of 2 x 10 6 x 0.035 ~ 70000 
unweighted events would have been produced. 

Moreover the following set of cuts has been applied: 

M u - d > lOGeV, | cos#(u(J),beam)| < 0.9, cos 6{u,d) < 0.9, E u{d) > 20GeV. 
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Both programs are equally fast and the run of 2 x 10 6 MC points costs no more 
than a few CPU minutes on DXPLUS@cern.ch. 



e e + 



In order to demonstrate the ability of PHEGAS/HELAC to deal with more complicated 
processes we give here results on four-fermion plus a gamma production. The results 
compare very well with the results presented in reference [[HJ form WRAP and 



RACOONWW [pT|] as is is shown in the following table: 



e~e + — > 


WRAP 


RACOONWW 


PHEGAS/HELAC 


udpTv^ 


75.732(22) 


75.647(44) 


75.683(66) 




78.249(43) 


78.224(47) 


78.186(76) 




28.263(9) 


28.266(17) 


28.296(22) 




29.304(19) 


29.276(17) 


29.309(25) 


udscq 


199.63(10) 


199.60(11) 


199.75(16) 



We refer to reference [10] for details on parameters and cuts used for this compu- 
tation, as well as an extensive comparison among the three generators based on 
differential distributions, which shows a very good technical agreement. 

gg -f bbbbW~W + 

The reason we have chosen such a process is twofold: in first place this is a challeng- 
ing process, from a computational point of view, and secondly this is a nice example 
to demonstrate the ability of PHEGAS/HELAC to deal with QCD processes. Moreover 
its study is important as a background of tiH production |12| . The results of the 
computation are summarized as follows: 



MC points 
w > 


result 
(fb) 


error 
(fb) 


efficiency 

(%) ' 


efficiency 

w > (%) 


99442 


4.716 


0.024 


3.3 


33 



The results refer to an energy yfs = 500 GeV and to 1 x 10 6 MC points. To give an 
idea of the complexity of the computation, the number of Feynman graphs for this 
process is 960, without taking into account electroweak contributions from Z and 7 
intermediate states. Parameters used are gQco = 1, mtop = 175 GeV and T top = 1.5 
GeV. Moreover the following set of cuts has been applied: 



M q „, > 20GeV, E q > 20GeV, | cos 8(q, beam)] < 0.9, 
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Figure 2: Differential distribution of the invariant masses mt 1 _w+ and mi_ w - for the 
process gg — > bbbbW + W~ . 



where q, q' refer to any quark or anti-quark of the final state. Finally in Fig|| we 
show the distribution of the invariant masses of b — W + and b — W~ pairs, exhibiting 
the expected peak at m top along with the non-resonant QCD corrections. 

In conclusion PHEGAS/HELAC provides an automatic and efficient computational frame- 
work to perform cross section evaluation and event generation for arbitrary scattering 
processes. 



Appendix 

In this appendix we describe the limits on t and Q\. The expression for the t invariant 
is given by 

t=(Qi- q 2 ) 2 = mj + Ql- ^(Q 2 + Ql- Ql) + ^-\q 2 \ cosfl 

with 

A = A(Q 2 , Ql, Ql) =Q A + Q\ + Q\- 2Q\Q 2 - 2Q\Q 2 - 2Q\Q\ 
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The limits for t can be found by maximizing (minimizing) t± given by 
t ± = ml + Ql- ^(Q 2 + Q\- Ql) ± ^\q 2 \ 

In order to find the maximum of t+ we study the function dt+ / dQ\ in the region Q\ min < 
Qi<(Q- Q2) 2 . Since 

dH + - -40 2 « 2 A- 3 / 2 M < 

d(Qir ~ q 

and 

we just consider two cases (\q 2 \ 7^ 0): 

1. 9t + /9Q 2 | Q 2 =Q 2 mm < in which case t max = t + , max = t+(Q\ = Q\, min ), and 

2. dt + / 8QI\q2 = q2 > in which case one can easily derive t max = t + {Q\ = x_) with 



Following the same reasoning we find for the lower limit on t that 

1. dt-/dQl\ Q 2 =Q 2 min > in which case t min = t- min = t-(Q\ = Qi >min ), and 

2. dt^/dQ\\Q2 = q2 < in which case one can easily derive t min = t-{Q\ = x+) with 

x + = Q 2 + Ql + 2QQ 2 1 ~ E * /Q 

V« 

The limits for the Q 2 -integration for given t can now be fixed by the condition | cos#| < 1 
or equivalently 

n(Q 2 ) < 

with 

n(Q?) = (t-Q\-m 2 2 + ^-(Q 2 + Qi- Ql)) 2 " (^) 2 A 

If Hi < Vi are the two roots of the polynomial TI(Qf ) then we have 

I. For a > , y_ < Q\ < y+, with y_ = max(yi, Q 2 jmin ) and y + = mm(y 2 , Q\ max ) 



2. For a < we have to satisfy two conditions Q\ < yi or y 2 < Ql and Q\ min < 
2 ^ /n2 
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